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ABSTRACT 

The factors affecting vortex growth in convectively stable protoplanetary 
disks are explored using numerical simulations of a two-dimensional anelastic- 
gas model which includes baroclinic vorticity production and radiative cooling. 
The baroclinic feedback, where anomalous temperature gradients produce vor- 
ticity through the baroclinic term and vortices then reinforce these temperature 
gradients, is found to be an important process in the rate of growth of vortices in 
the disk. Factors which strengthen the baroclinic feedback include fast radiative 
cooling, high thermal diffusion, and large radial temperature gradients in the 
background temperature. When the baroclinic feedback is sufficiently strong, 
anticyclonic vortices form from initial random perturbations and maintain their 
strength for the duration of the simulation, for over 600 orbital periods. 

Based on both simulations and a simple vortex model, we find that the local 
angular momentum transport due to a single vortex may be inward or outward, 
depending its orientation. The global angular momentum transport is highly 
variable in time, and is sometimes negative and sometimes positive. This result 
is for an anelastic gas model, and does not include shocks that could affect angular 
momentum transport in a compressible-gas disk. 

Subject headings: accretion, accretion disks, circumstellar matter, hydrodynam- 
ics, instabilities, methods: numerical, turbulence, solar system: formation 
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1. Introduction 

The baroclinic term is a source of vorticity in tlie vorticity equation, and is derived by 
taking the curl of the pressure gradient in the Navier-Stokes equation, 

V X (-^Vj>^ = ^Vp X Vp, (1) 

where p is the pressure and p is the density. The baroclinic term is nonzero when pressure 
and density gradients are not aligned. 

An intuitive example of baroclinicity is the land-sea breeze, which is initiated when air 
temperatures above the land rise more than over the nearby ocean. The warm air over the 
land expands, isobars rise relative to those over the ocean, and consequently the isobars tilt 
towards the ocean. At the same time, the colder air over the ocean has a higher density than 
over the land, so the isopycnals tilt towards the land. The tilting of isobars and isopycnals 
in opposite directions is a baroclinic source of vorticity, which causes a circulation in the 
vertical plane that blows from the ocean to the land near the surface. Thus the potential 
energy of the tilted isopycnals is converted into kinetic energy of the land-sea breeze, which 
dissipates through surface fricti on and reduc es the land-sea temperature contrast through 



temperature advection (see, e.g. lHoltonll2004l ) 



A related concept is the baroclinic instability, which is of central importance to the 
production of vortices and Rossby waves at midlatitudes. Here the decrease in solar insolation 
from equator to pole causes colder temperatures, and consequently higher density, at the 
surface at higher latitudes; thus the isopycnal surfaces are tilted towards the equator. A 
system with tilted isopycnals has more potential energy than one with level isopycnals, just 
like an inclined free surface has more potential energy than a level one. This potential 
energy is available to processes which can flatten out the isopycnals. For example, vortices 
in the atmosphere and ocean convert the potential energy of the inclined isopycnals to the 
kinetic energy of their meso-scale motion. Vortices flatten the isopycnals by transferring 
heat poleward through their mixing action. 

The baroclinic instability is so-named because of the tilted isopycnals, but the physics 
is fundamentally different from the land-sea breeze: in the land-sea breeze the circulation is 
in the vertical plane and caused directly by the baroclinic term, i.e. by non-aligned density 
and pressure gradients in the vertical; in the baroclinic instability the isopycnals are titled in 
the vertical, but the vortices are in the horizontal plane, so they could not be produced by 
the baroclinic term directly. Rather, the tilted isopycnals present an unstable configuration 
which is ripe for processes which can convert the potential energy to kinetic energy, much 
like how an avalanche levels out a steep incline of snow. 
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The baroclinic processes discussed in this paper for a protoplanetary disk are similar 
to the land-sea breeze, but in radial geometry. Due to the gravity and radiation of the 
central star, the density, temperature, and pressure of the disk's gas all decrease radially. 
Any azimuthal variations in temperature (and thus density or pressure by the ideal gas law) 
would lead to an increase in vertical vorticity due to the baroclinic term ([1]). The focus of this 
work is the baroclinic feedback, where a vortex enhances azimuthal temperature gradients 
to reinforce the vortex itself. Under the right conditions, the baroclinic feedback strengthens 



planetary formation, as they are efficient at collecting particles from the disk ( 


Tanea et al. 


1996; 


Johansen et al. 


2004; 


Baree & Sommeria 


1995; 


Klahr & Bodenheimer 


20061). The high 



density of solids in the vortex would speed the formation by core accretion, which is so slow 
i n the rest of th e disk that it may not be a feasible theory of planetary formation there 
( IWetherilll Il990l ). A vortex which collects solids is also a potential site of gravitational 
instability, where the matt er is dense enough that it simply collapses into a planet through 
gravitation self- attraction (IBossI 119971 ). 



Our study was motivated by lKlahr fc Bodenheimer! ((20031), who investigated the effects 
of baroclinicity in a radially stratified disk using a finite difference model of the compressible 
Navier-Stokes equation combined with a radiative transfer model. They found the baroclinic 
instability to be a source of vigorous turbulence which leads to the formation of long-lasting 
vortices and positive angular momentum transport. Barotropic simulations, where the en- 
tropy (temperature) is constant in the radial directio n did not deve lop turbulence, even with 
large initial perturbations. To explain these results, iKlahr performed a local linear 

analysis for a disk with constant surface density, and found that modes do not grow if the 
growth time of the instability is longer than the shear time. 

The issue of whether the baroclinic instability is a mechanism for nonlinear growth and 
the formation of vortices has been a recent source of debate. Johnson and Gammie are critical 



of the findings of lKlahr fc Bodenheimerl (120031 ) and lKlahrl (120041 ). Their linear analysis found 
no ex ponentially growing instabil ities, except for convective instabilities in the absence of 
shear (j Johnson fc Gammiell2005al ). Furthermore, they use a shearing sheet numerical model 
to show that disks with a nearly-Keplerian rotation profile and radial g radients on the order 
of th e disk radius are stable to local nonaxisymmetric disturbances (jjohnson fc Gammie 
2006h . 



The goal of this study is to understand the effects of baroclinic instabilities and ra- 
diative cooling on the generation of turbulence, vortex formation, and vortex longevity in 
protoplanetary disks. One o f our motivations is to s hed li ght o n the conflicting observatio ns 
of baroclinic instabilities by lKlahr fc Bodenheimerl (120031 ) and lJohnson fc Gammid (120061 ). 
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This work is presented in two parts. Part I, which precedes this article, presents the 
equation set, details of the numerical model, and results of the small-domain simulations, 
which are used to study the process of vortex formation. This paper. Part II, explores the 
parameters which affect the baroclinic feedback during the growth phase of the vortices; 
these simulations use the larger quarter annulus domain and are run for hundreds of orbital 
periods to observe the long-term behavior of the vortices. We begin with a quick review of 
the equation set in §2] The results in ^ discuss the evolution of a typical simulation, the 
process of the baroclinic feedback, the Richardson number as a diagnostic, and the alpha 
viscosity. In ^ we discuss the angular momentum transport in our simulations, which is 
highly variable and depends on the orientation of individual vortices. In [|5]we conclude that 
the baroclinic instability is an important mechanism for vortex generation and persistence, 
and review the conditions which affect the instability. For conciseness there is little repetition 
between Parts I and II, so the reader is advised to read both together. 



2. Description of the Equation Set 



The model equations are described fully in Part I of this work, and are only briefly 
reviewed here. They model an anelastic gas, which filter out pressure waves that restrict 
the timestep of the numerical m odel but do no t im pact the physics of interest here. Our 
equations set is similar to those in lBannoru (119961 ) and lScinocca fc Shepherd! (119921 ). which are 
anelastic models of the atmosphere derived from conservation of momentum, conservation 
of mass, the second law of thermodynamics, and the ideal gas law. Our equations use 
two-dimensional polar coordinates (r, 0) where temperature and density are stratified in the 
radial direction. Variables such as the vertical component of vorticity (, streamfunction \E^, 
potential temperature 6, thermal temperature T, surface density S, and Exner pressure vr 
are written as the sum of a background and perturbation term, e.g. 6 = do{r) + 6'{r,(f),t), 
where the background functions only vary radially. 



The model equations are 



1 d f r d^' 



'dt 



r dr \ Eq dr 

Cp Ottq 89' 
r dr d(j) 
6' 



1 d^^' 

2/-I 



T 



(2) 
(3) 
(4) 



The first is the relationship between the perturbation streamfunction $' and the perturbation 
vorticity C'; the second and third are prognostic equations for perturbation vorticity C,' and 
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perturbation potential temperature 9' . Radial and azimutlial velocities u = (u, v) are related 
to the streamfunction by Squ = —V x \l/z. Other variables include the radiative cooling time 
r, specific heat at constant pressure Cp, time t, viscosity z/g, thermal dissipation Ke, vertical 
unit vector z, and the Jacobian d{a, h) = [dr-ad^h — d^adr-h)/r. 



The baroclinic term 



Cpd7fod9^ 
r dr d(j)^ 

is a central focus of this paper. It is the only source term in the vorticity equation ([3]), and 
it plays an important role in the baroclinic instability, as one might expect. Most people 
are familiar with the baroclinic term using density and pressure, shown in eqn. ([1]). This 
operation in terms of our variables is 

^ d f ^, diTo \ Cp 39' diTQ 



^ ^ r d(j) \ ^ dr ) r d(j) dr ^ ^ 

The factor driio indicates that the baroclinic feedback should be strengthened if (9,.7ro is large, 
i.e., if radial pressure gradients are large. But 

dr dr dr dr ' 

so we expect large radial temperature or density gradients to strengthen the baroclinic feed- 
back. 

The radiative cooling term 9'/ t diffuses the perturbation potential temperature equally 
at all scales with an e-folding time of r. The two Laplacian terms, UeV'^C and Ke'V'^9', 
diffuse potential vorticity and potential temperature at fastest at the highest wavenumbers. 
They were added to the numerical model to dissipate energy for numerical stability. In 
the nondimensionalized version of the equation set the u coefficients are replaced with the 
Reynolds and Peclet numbers, 

7-2 7-2 

Re = Pe 



where the length scale Lsc and time scale tgc are described in Part I. 



3. Results 



The simulations discussed in this paper vary parameters such as the radiative cooling 
rate, the background temperature and surface density gradients, and the Peclet number 
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(Table [T]). These simulations capture the salient features of the physics of the anelastic 
equation set. The topic of Part 1 was vortex formation, and thus used a smaller domain 
for only five orbital periods. Here we are interested in vortex growth and longevity due 
to the baroclinic feedback, and have chosen a larger domain and durations of 300 to 600 
orbital periods. (This is 6,200 to 12,400 years for a solar-mass star.) The simulations were 
performed on the quarter annulus with a radial extent from 5AU to lOAU and a resolution 
of 256 X 256 and 512 x 512 gridpoints. 

The background surface density and background temperature are constant in time and 
are power functions in the radial. 



So(r) = a — 



To(r) =ci — 



(8) 



where rj„ = 5AU is the inner radius of the annulus. The coefficients are a = lOOOg cm~ , 
c = 600K for the quarter annulus domain; b and d are varied and shown in Table [TJ For 
example, for simulation Al, the background surface density varies from lOOOg cm"^ to 350g 
cm~^ and the background temperature decreases radially from 600K to 150K. This range of 
temper atures can only be achieved in a realistic disk when the radius ranges from 1 AU to 
10 AU flBosslll998l ). We have artificially enhanced the radial temperature gradient in order 
to compensate for the lower resolution of our global simulations. We have demonstrated in 
paper I using a higher-resolution local simulation that more realistic temperature gradients 
can still produce vortices. Most simulations were run to 300 orbital periods, measured as a 
full (27r) orbit at Tmid = 7.5AU. This is 6,200 years for a solar-mass star. 



Thermal temperature T and potential temperature 6 are related by 



e = T 



Po[rir. 
P 



R/cp 



T 



vr 



(9) 



where R is the gas constant and tt is the Exner pressure. All results in this paper are 
expressed in terms of the thermal temperature T in order to compare to observations. The 
potential temperature is a measure of entropy. If entropy increases r adially (d9n/dr > ) 
then the disk is convectively stable — this is the Schwarzschild criterion (jSchwarzschildlll958l ). 
If the entropy gradient as accompan ied by differential rotation, the Solberg-H0iland criterion 
( iTassoull I2OOOI : iRiidiger et al.l |2002| ) is used to test convective stability (see Results section 
of Part I). For the simulations presented in this paper, the Solberg-H0iland value is positive 
(0.035-0.299 years"^), indicating that they are convectively stable. 

The initial condition for the perturbation temperature is shown in Fig. [1] It is created 
with a specified wavenumber distribution in Fourier space, transformed to Cartesian coordi- 
nates, and interpolated to the Fourier- Chebyshev annular grid (see Part I, section 3.2). The 
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initial vorticity perturbation is created in a similar fashion. The magnitude of the initial 
conditions is 25% of the maximum of the background function. 

The small domain simulations in Part I (r G [9.5, 10], G [0,7r/32]) required a much 
smaller initial perturbation to initiate vortices — a temperature perturbation of only 5%, 
and an initial vorticity perturbation of zero. This is possible because the small domain 
is a higher resolution relative to the background shear. The sensitivity analysis in Part I 
showed that smaller initial perturbations are required to initiate vortices with progressively 
higher resolution and Reynolds number. The same is true of the background temperature; 
the quarter annulus domain uses higher temperatures (c = 600K) and steeper gradients 
[d = —2) than the small domain simulations. Again, the sensitivity analysis in Part I 
showed that at higher resolutions vortices can be formed with progressively cooler disks and 
shallower background gradients. 

The evolution of a typical simulation can be described as follows: The initial distribution 
of vorticity shears due to the differential rotation of the nearly-Keplerian rotational profile 
(Fig. [2]). Even at these early times, the perturbation vorticity and perturbation kinetic 
energy grow due to the baroclinic term (FigE]). After about five orbital periods, anticyclonic 
vortices begin to form, and by ten orbital periods the domain is populated by numerous small 
anti-cyclones. Cyclonic (anti-cyclonic) fluid rotates in the same (opposite) direction as the 
background fluid, and is denoted by positive (negative) vorticity perturbation in the figures. 
It is well-known that anti-cyclones can be long-lived in a Kepler ian disk, while cyclones shear 



out into thin filarn ents that eventually dissipate away fiGodon fc Livid Il999l : lMarcuslll990 



Marcus et al.ll200d ). An anticyclonic vortex has a positive azimuthal velocity at small inner 
radii and a negative azimuthal velocity at large outer radii. This means that anticyclonic 
vortices can smoothly match the background shear flow, and therefore extract energy from 
the Keplerian shear. Cyclonic vortices cannot smoothly match the background shear flow 
and are therefore sheared apart. 

After the initial period of vortex formation, the vortices merge and grow in strength 
(Figs. [3], Hj). This merging behavior is similar to the merging of like-signed vortices in two- 
dimensional isotropic turbulence, which transfers energy from smaller to larger scales (the 
inverse cascade). However, in shearing flows vortices do not merge as readily and must be 
sufficiently close in the radial direction. It is not at all clear that this merging of vortices can 
occur in a fully three-dimensional disk if the initial radial vortex scale is small compared to 
the disk scale height. On the other hand, if vortic es primarily form ori the up per and lower 



surfaces of a vertical stratified disk as found by iBarranco fc MarcusI (120051 ) . then it may 



be possible for small-scale vortices to merge in these surface layers. Further discussi o n and 



images of vortex merger, longevity, and distribution can be found in iGodon fc Livid (119991 ) 
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and lUmurhan fc RegevI (120041 ) . 

There is a clear "sandwich" pattern of temperature perturbations around each vortex 
(Fig. Hj): the vortex advects warmer fluid towards the outside of the disk and cooler fluid 
towards the inside of the disk. In the sandwich analogy, the temperature perturbations are 
the bread and the vortex is the meat between the bread. These perturbations have azimuthal 
temperature gradients that play a role in the baroclinic feedback. 



3.1. Baroclinic vorticity production 

The model equations for vorticity and temperature perturbations are coupled by the 
baroclinic term in the vorticity equation ([3]) and the advection term in the temperature 
equation (j4]). This coupling is required to support long-lived vortices; without it, vorticity 
and temperature perturbations simply decay to zero. 

The baroclinic feedback operates as follows: 

1. Azimuthal gradients in the perturbation temperature field, dQ' jd(^^ make the baroclinic 
term in the vorticity equation non-zero. 

2. The baroclinic term is a source of vorticity which strengthens anticyclonic vortices. 

3. Vortices stir the fluid, moving warm fluid from the inner annulus outward and cool 
fluid from the outer annulus inward. 

4. This local advective heat transport enhances azimuthal temperature gradients, 89' /dcj), 
completing the feedback cycle. 

In order to show that the vortex growth is indeed due to this baroclinic feedback, the 
baroclinic term was turned off at various times in simulation set B (Fig. [5]). In all of these 
trials, perturbation vorticity and kinetic energy drop off immediately when the baroclinic 
term is turned off. This is particularly striking at t = 10 and t = 100, when vortex strength 
is growing quickly in the reference simulation. The kinetic energy in these plots is computed 
from the perturbation velocity fields. 

The rate of thermal dissipation, r, plays a crucial role in the formation and growth 
of vortices. Fig. E] shows that there are two distinct stages in these simulations: vortex 
formation, from t = to about t = 5 orbital periods, and vortex growth, which occurs after 
t = 5. During vortex formation, small thermal dissipation (large r) allows the strongest 
vortices to form. That is because the initial temperature perturbation dies off quickly when 
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thermal dissipation is large, so that azimuthal temperature gradients gradients are smaller 
and the baroclinic term produces less vorticity. This is not yet the baroclinic feedback 
because steps 3 and 4 are missing — it is just baroclinic vorticity production from the initial 
temperature gradients, which are steps 1 and 2. 

Once vortices form they advect fluid about them (step 3), creating the distinctive "sand- 
wich" pattern of cool (warm) temperature perturbations to the inside (outside) of the vortex, 
as shown in Fig. HI These temperature perturbations create local azimuthal temperature 
gradients (step 4), completing the cycle of the baroclinic feedback. Sometime after five 
orbital periods, the vortices have formed and the simulation transitions from the vortex 
formation stage to the vortex growth stage. Now that the baroclinic feedback is operating, 
thermal dissipation has the opposite effect than at early times (Fig. E]). If the disk cools 
quickly (r small), then the warm and cool temperature perturbations can remain tight about 
each vortex, so that 89' /d(j) in the baroclinic term is large, and the baroclinic feedback is 
strong. If the disk cools slowly (r large), the perturbation temperature responds sluggishly 
to mixing by vortices, 89' /del) is small, the baroclinic feedback is weak, and vortices simply 
dissipate away (Fig. [7j). Quantitative measures of disk activity like kinetic energy, and max- 
imum temperature and vorticity clearly show that the strength of the feedback and rate of 
growth of vortices is strongly dependent on r (Fig. [6]). In simulations where the radiative 
cooling rate was sufficiently fast, the baroclinic feedback counters dissipation and vortices 
remain strong and coherent for hundreds of orbital periods (Fig. Hj). The longest running 
simulation, T3 where d = —0.75, ended at 600 orbital periods, at which point all of the 
vortices had merged into a single anticyclonic vortex. 

There are two dissipative terms in the temperature equation (jlj): the Laplacian term 
KeV^6'', which dissipates most quickly at small scales; and the radiative cooling term —9'/t 
which dissipates equally at all scales. Can the Laplacian term play the same role as the radia- 
tive cooling term in the baroclinic feedback? Simulations Pel-Pe3, where Pe ranges from 10^ 
to 2 X 10^, show that the Laplacian term can indeed play that role (Fig. [8]); higher thermal 
diffusion (smaller Peclet number) produces a stronger baroclinic feedback. Higher diffusion 
produces warm and cool areas around each vortex which are more localized azimuthally, and 
therefore have larger azimuthal temperature gradients (Fig. The azimuthal temperature 
gradients then produce more vorticity through the baroclinic term (step 1 of the baroclinic 
feedback) . 

Other simulations explore the role of background temperature (T1-T5) and background 
surface density (D1-D3). Larger background temperature gradients in simulations T1-T5 
result in larger and stronger vortices (Fig. [TOl) . Quantitative measures such as the kinetic 
energy, maximum vorticity, and maximum temperature all grow faster with larger tempera- 
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ture gradients (Fig. [TT]) . The evolution of these quantities does not change as the background 
density gradient is varied (simulations D1-D3, Fig [T2|) . It is somewhat surprising that the 
baroclinic feedback responds strongly to the background temperature gradient but not the 
background density gradient, when these gradients seem to be on equal footing in the baro- 
clinic term (see eqn. [7]). The background temperature gradient is a source of available 
potential energy that can be transformed into the kinetic energy of vortex motion as the 
vortices transport heat from the hot inner disk to the cold outer disk. This nonlinear heat 
advection cannot be captured in a linear stability analysis. Since the surface density is time- 
independent in our anelastic model, the background surface density gradient cannot provide 
a source of potential energy for vortex formation. 



3.2. Richardson number 



Several previous st udies have used the Richardson n umber to characterize instabilities 
in protoplaneary disks (jjohnson fc Gammid l2005al . l2006l ) , and we compute the Richardso n 
number here for cora parison. We believe that the Solberg-H0iland criterion (ITassoull (120001 ): 
Riidiger et al.l (|2002| ). also see §4 of Part I), which was specifically created for differentially 
rotating astrophysical fluids, is the best way to judge whether a disk is convectively unstable. 
For the simulations presented in this paper, the Solberg-H0iland values are positive (0.035- 
0.299 years"^), indicating that they are all convectively stable. However, the Richardson 
number also provides useful information about the instability. We found that the baro- 
clinic feedback is stronger (i.e. vortex growth is faster) in simulations with more negative 
Richardson numbers. 

The Richardson number is often evoked in geophysical turbulence to quantify the re- 
lationship between stratification and shear. For the atmosphere this dimensionless ratio is 
typically 



Ri( 





_g_dp_ 


N\z) 


p dz 




fduV 







(10) 



where N{z) is the local Brunt- Vaisaa buoyancy frequency, u is the hori zontal velocit y, z is 
the vertical coordinate, p is the density, and g is the gravitational force (jTurnerlll973l ). The 
numerator A^^ gives the strength of stratification, where iV^ is negative for a convectively 
unstable fluid, is positive and small for weakly stable stratification, and is positive and large 
for strongly stable stratification. The denominator gives the strength of the shear. 
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In our equation set the Richardson number is 



= -T^ 2 = ^H^- (11) 



dr j \ dr 

By comparing the Richardson number (Fig. [T3!l with kinetic energy or maximum vorticity 
(Fig. [TTl) for simulations T1-T5, it is clear that the Richardson number is an excellent way 
to predict the strength of the baroclinic feedback. When Ri > (T4 and T5, where d = —0.5 
and —0.25), kinetic energy and vorticity simply decay away. When Ri < (T1-T3, where 
d = —2 — 0.75), the baroclinic feedback operates and kinetic energy and vorticity grow. In 
fact, the simulation with the most negative Richardson number (Tl, where d = —2) also has 
the fastest vortex growth. 



Johnson fc Gammid (120061 ) found that disks with a nearly-Keplerian rotation profile and 
radial gradients on the order of the disk radius have Ri > —0.01, and are stable to local 
nonaxisymmetric disturbances. Our simulations are not restricted to this Ri > —0.01 crite- 
rion, as simulations T1-T3 have quickly-growing instabilities but have Richardson numbers 
in the range of —5 x 10"^ to —5 x 10~^. The most likely difference between the two models 
that accounts for this disagreement is that our simulation allows small initial temperature 
perturbations to evolve into strong local vorticity perturbations that can produce stable vor- 
tices. This initial evolution can only occur if the viscous dissipation is sufficiently low (high 
Reynolds number). 



3.3. Alpha viscosity 



Protoplanetary disks are often described by the dimensionless number a, which is used 
to parameterize an effective viscosity u = aCgHp where Hp is the vertical pressure scale height 
of the disk and Cg is the local sound speed. This simple description was used to calculate the 
density structure, temperature structure, and mean components of laminar and turbulent 
gas fl ow in a disk (IShakura fc Sunyaevlll973l . iLynden-Bell fc Pringlelll974l . iLin fc Papaloizou 
1980h . 



Alpha viscosity, rather than Reynolds number, is commonly reported in the astrophysi- 
cal literature to characterize the dissipation of energy in the disk. If the pressure scale height 
is scaled as Hp = Cs/flo, where Qq is the background angular velocity, the alpha viscosity 
can be calculated as a{r) = i'eQo{r)/c1. In our anelastic model, this measure cannot be used 
directly because Cg » \u'\ and pressure waves are temporally constrained to adjust instan- 
taneously. In order to compare the alpha viscosity with other protoplanetary disk models 
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we report the ratio of the alpha viscosity to the Mach number squared, 



a. 



u 



Relu 



(12) 



where Ma = lu'l/c, and tildes indicate nondimensionalized variables. 



Azimuthal averages of a/Ma^ for all simulations are between 10 ^ and 10 ^ (Fig. [T^ . 
Mach numbers of 0.0 1 or 0.1 would produce corre sponding alpha viscosity ranges of 10"^- 



10 



a 



^9 and lO^MO-^. 

= 10"^ to 10""^ for their two-dimensi onal simulations with radial temperature gradients, 



Klahr fc Bodenheimerl (120031 ) report Mach numbers of 0.03 to 0.3 and 



which have resolutions of 62^ and 128^. I Go don &: Livid (119991 ) report a viscosity parameter 
a = 10~^ and 10~^ for their 128^ and 256^ simulations, respectively. In general, higher 
Reynolds numbers (and thus smaller alpha viscosity) can be achieved with higher resolution. 
Our simulations have slightly higher resolution (256^ and 512^) than other studies, and the 
effective Reynolds number is also higher (~ 10^) due to the use of hyperviscosity (see Part 
I, section 3). These characteristics contribute to a lower effective alpha viscosity, so that 
our results include numerous fine, small-scale structures such as layers of filaments swirling 
around the vortices. 



4. Angular Momentum Transport 

The transport of angular momentum is of critical interest in the study of protoplanetary 
disks. The traditional view of disk evolution is that angular momentum is transported 
outward as mass is transported inward. In Keplerian motion gas near the star has a faster 
angular velocity than the gas further out. Turbulence in the gas creates an effective viscosity, 
so that faster moving gas in the inner disk will speed up slower gas in the outer disk, and the 
outer fast gas will tend to slow down the inner gas. Thus angular momentum is transported 
outward. As the inner gas slows down, it is no longer rotationally supported at that orbit, 
and falls toward the star to gain speed. Thus mass is transported inward. Similar arguments 
can be made for particle collisions, which would enhance this process. 

The theory of outward angular momentum transport is based on azimuthally uniform 
dynamics in a viscous disk. Turbulence and coherent structur es may have radically d i fferent 



effects, and is currently a topic of intense scientific interest. iKlahr Sz Bodenheimerl (120031 ) 



report that, just like in laminar fiow, turbulence in baroclinic disks tr ansports angula r mo 



mentum outward and mass inward, while releasing potential energy. iLi et al.l (12001! ) used 
a finite volume model of the compressible Euler equations to model Rossby waves and vor- 
tex generation, and found that individual vortices transport angular momentum outward. 
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Johnson fc Gammid (l2005bl ) also found positive angular momentum flux in their compress- 
ible shearing sheet model when they used strong initial vorticity perturbations to trigger 
vortex formation. 

As a simple example, consider locally Cartesian coordinates in the radial and azimuthal 
direction. A slanted vortex of this local coordinate system could have the stream function 



^' = aexp\-\ f +i3r<j)+{-) . (13) 




Each streamline is a rotated ellipse centered at the origin with radial extent tq and azimuthal 
extent (po. The angle of the ellipse is only affected by jS. This vortex is superimposed on 
some background flow, so the perturbation velocities in locally Cartesian coordinates are 
u' = —d(f,'$' and v' = dr"^'. The angular momentum transport of this vortex is 



nVd0 = -/3aVoSo^exp U (^^^ (-4 + /J^gro^) j . (14) 



The sign of this quantity depends only on f], the angle of the vortex. Positive and negative 
vortices with the same orientation have the same angular momentum transport, as the sign 
of the vortex only affects a, a squared quantity in f|T^ . This indicates that it is only the 
orientation of the vortex within the flow that affects whether momentum travels towards the 
inside or outside of the disk; the direction of rotation of the vortex is inconsequential. 

Our simple analytic example is shown in Fig. [15] for j3 = —0.5 (top left) and (3 = 0.5 
(top right), where the other constants are 0o = Ij '''o = 2, and a = 1. The bottom panels 
show vortices with similar orientations in the full numerical model. Clearly, the direction 
of the angular momentum transport only depends on the angle of the vortex, as in the 
analytic example. These vortices are not from the simulations in Table [H but are from short 
simulations that were speciflcally designed to produce these orientations. 

What is the effect of vortices on angular momentum transport when they are imbedded 
in a turbulent flow populated with filaments and other interacting vortices? To investigate 
this, the angular momentum transport, Hqu'v', was recorded using azimuthal and global 
averages in the numerical model. In typical simulations, like Al, the angular momentum 
transport is highly variable in space and time (Fig. [T6l) . Specifically, the global angular 
momentum transport cycles chaotically between positive and negative periods as the vor- 
ticity field evolves. The angular momentum transport in these simulations is infiuenced by 
the interaction of numerous vortices and vorticity filaments, which is much more compli- 
cated than the single vortex case. We would expect that individual vortices within this fiow 
would contribute angular momentum based on their orientation, and that these individual 
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contributions could be summed to find the angular momentum transport. However, a sepa- 
rate study of a small number of vortices and filaments in the flow would be required to say 
conclusively. 



Based on the simulation results in Fig. [161 we conclude that the total angular momen- 
tum transport in an anelastic-gas turbulent disk with vortices and vortex filaments may be 
inward or outward, and can vary locally in the disk depending on the orientation of the vor- 
tices. In contrast, studies of compressible- gas disks have all found that vortices transport an- 
gular m omentum outwards (IKlahr &: Bodenheimerll2003l : iLi et al.ll200ll : iJohnson &: Gammie 
2005bl ). Compressible-gas models include acoustic waves, which are filtered out of our anelas- 
tic model. Shocks produced by acoustic waves in these studies may orient the vortices uni- 
formly so that angular momentum is transported outwards, or transport angular momentum 
by other means. 



Conclusions 



In this study we are interested in exploring the necessary conditions for vortex formation 
in an anelastic protoplanetary disk model that includes baroclinicity and radiative cooling. 
We have shown that long-lived vortices can be formed by initial random temperature per- 
turbations through the mechanism of the baroclinic instability. Vortex production must 



compete with the strong inhibiting effects of Ke 



jlerian shear, an effect observed by other 



authors (IBracco et al.lll999l . iGodon &: Livid Il999l ). Only anticyclones survive in Keplerian 



disks, while cyclones shear out and diffuse awa y. Many previous studies do not include baro 



clinic effects due to a lack o f thermodynamics (IBracco et al.lll99g 



Umurhan 



fc RegevI 



2004 



Johnson &: Gammie 2005b ) or an assumed polytropic relation ( Godon &: Liviol 1999 ). and 



therefore only model decaying turbulence from an initial vorticity distribution. 

In the baroclinic feedback, local azimuthal temperature gradients produce vorticity 
through the baroclinic term in the vorticity equation. This strengthens vortices, which 
advect the surrounding thermally stratified gas, producing stronger local temperature gra- 
dients. 

The baroclinic feedback can only operate once vortices have formed, as a coherent 
vortex is required to produce the "sandwich" pattern of warm and cold gas about each 
vortex. In our simulations two distinct stages can be seen: vortex formation, where the 
initial temperature perturbation rapidly decays; and vortex growth, where the baroclinic 
feedback takes effect and both perturbation vorticity and perturbation temperature grow for 
the rest of the simulation. 
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The conditions required for the barochnic feedback are: (1) a sufficiently large radial 
temperature gradient in the background stratification; and either (2) a fast radiative cool- 
ing time; or (3) high thermal dissipation (i.e. small Peclet number). If the background 
radial temperature gradient (condition 1) is too small, advection by the vortices does not 
strengthen local azimuthal temperature gradients. Both conditions 2 and 3 allow tempera- 
ture perturbations to track vortices so that the structure of the vorticity and temperature 
fields are strongly coupled. The difference between the two mechanisms is that thermal 
dissipation smoothes out small scale features, resulting in large-scale thermal perturbations, 
while radiative cooling affects all scales equally and produces smaller thermal perturbations. 
Varying the background surface density gradient had no effect on the strength of the baro- 
chnic feedback. 



One of the goals of this study was to see if the baroclinic instabilities found by lKlahr fc Bodenheimer 



( I2OO3I ) can be reproduced in an anelastic equation set with simplified dynamics. They found 
that if the background radial entropy gradient is zero — this turns off all baroclinic effects — 
then the initial vorticity perturbation just decays away (their Model 2). When entropy varies 
radially so that the temperature T ~ r~^, the flow becomes turbulent within a few orbits 
and vortices are formed (their Models 3-6). Our results (condition 1, above) agree with 
this result, and furthermore show that vorticity grows faster with steeper background radial 
temperature profiles. 

Our conditions 2 and 3 state that thermal dissipation is required for the baroclinic 
instability. Indeed, in our simulations where both forms of thermal dissipation were suffi- 
ciently slow, the vorticity dies of f after the initial vortex forrn ation. This requirement is in 



disagreement with the findings of iKlahr fc Bodenheimerl (120031 ). as they "got rid of radiation 



transport" (i.e., there is no radiative cooling) for their 2D simulations. Because we use a 
simplified one-parameter radiative cooling model, we can see that t he baroclinic feedback 



strong ly depends on the cooling time r. It is not clear to us why iKlahr &: Bodenheimer 



(I2OO3I ) see vortex growth when radiative cooling is missing. Either there is some implicit 
thermal diffusion in their code, or their vortices are growing through a different mechanism 
than ours. 



The range of Richardson numbers at which we form vortices differs from lJohnson fc Gammie 



(120061 ). who find that simulations with Ri > —0.01 are stable to local nonaxisymmetric dis- 
turbances. Our simulations with —0.01 < Ri < form turbulent instabilities and vortices 
quite easily. A likely explanation for this difference is that our simulations have the large 
Reynolds numbers required to permit small initial temperature perturbations to evolve into 
strong local vorticity perturbations before they are viscously damped. 

The results of a model must be understood within the asymptotic regime where it is 
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valid. Our model assumes that the disk is thin and hydrostatically balanced in the vertical, 
so that only large-scale horizontal motions are consi dered. The vertical dyna mics, which 
we don't consider, can affect vortex stability as well: iKnobloch fc SpruitI (119861 ) argue that 
height variations must be included when discussing she ar instabilities in the disk, b ecause 
vertical gradients of azimuthal velocity are not small; iBarranco fc MarcusI (120051 ) found 
that columnar vortices are unstable to small perturbations, but that internal gravity waves 
naturally create robust off-midplane vortices (see discussion in conclusions of Part I). Indeed, 
the most serious limitation of this study is our assumption of a two-dimensional disk. If 
the initial baroclinic instabilities have radial scales that are small compared to the disk 
scale height, as is suggested by our local simulations presented in paper I, then the vertical 
stratification of the disk will likely play a major role in the nonlinear development and the 
longevity of vortices. Nevertheless, we believe that our two-dimensional simulations have 
served to identify physical processes that will likely play a role in a fully three-dimensional 
simulation. For example, if vortices primarily form on the upper and lover surfaces of disks, 
then the radiative cooling rates will likely be more rapid than if the vortices were buried in 
the optically thick midplane of the disk. We therefore expect that vortices confined to the 
surface of a disk could have longer lifetimes due to their ability to efficiently transport heat 
radially outwards. 

Our model is based on an anelastic mass conservation equation, which filters out acoustic 
waves. Thus shock waves and their potential interactions with vortices and angular momen- 



tum t ransport do not appear in our study, as they have in coin pressible gas models (ILi et al. 



200ll : iKlahr fc Bodenheimerll2003l : I Johnson fc Gammiell2005bl ). These studies all found that 
angular momentum is always transported outward, while we found that it may be inward or 
outward, and is highly variable in space and time. The obvious difference in the dynamics is 
the lack of shocks in our anelastic-gas model. This suggests that shocks play an important 
role in the transport of angular momentum in protoplanetary disks. 

The baroclinic feedback enhances vortices so that they can be long-lived. In our sim- 
ulations they survived for the duration of the longest numerical simulation — for over 600 
orbital periods (12,400 years) — and show no signs of decaying. This study shows that the 
baroclinic feedback is a viable mechanism for the generation and persistence of vortices in 
protoplanetary disks. In the baroclinic feedback, the background temperature gradient pro- 
vides a source of available potential energy that can drive the vortices indefinitely even in 
the presence of a finite rate of viscous dissipation. As vortices are efficient at collecting 
particles from the surrounding gas, they are a natural place for planets to form in the disk. 
If vortices are long-lived coherent structures in protoplanetary disks, as suggested by this 
work, they offer a way to overcome the difficulties presented by current planetary formation 
theories: the high particle concentrations in vortices speed up the core accretion process; 
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likewise, this high particle concentration could lead to gravitational instability. Both core 
accretion and gravitational instability are hindered in the majority of the disk where particle 
concentrations are low. Strong concentration of particles may require the vortices to grow 
large compared to the scale height of the disk so that they can extend through the midplane 
of the disk where most particles will reside. Fully three-dimensional simulations are therefore 
required to establish of relevance of vortices to planet formation. 
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Fig. 1. — Initial temperature perturbation T'. 
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pert Vort, t=0.5 orb per pert Vort, t=5 orb per pert Vort, t=10 orb per 




Fig. 2. — The perturbation vorticity (' in the quarter-annular computational domain for 
simulation Al. The time t refers to the orbital period in the middle of the annulus. At very 
early time (left), the vorticity is simply sheared by the background differential rotation. By 
five orbital periods (center) a few anti-cyclonic vortices begin to fold over, and by ten orbital 
periods (right) numerous small anticyclonic vortices have formed. 



- 22 - 



pert Vort, t=35 orb per pert Vort, t=51 orb per pert Vort, t=87 orb per 




pert Temp, t=35 orb per pert Temp, t=51 orb per pert Temp, t=87 orb per 




Fig. 3. — The perturbation vorticity (' (top) and perturbation temperature T' (bottom) for 
simulation Al, where the radiative coohng time is fast (r = 1). In this regime the barochnic 
feedback remains strong, and vortices remain strong for the full simulation. 
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pert Vort, t=300 orb per pert Temp, t=300 orb per 




V ■ 




Fig. 4. — Perturbation vorticity and temperature for simulation T3, where d = —0.75. De- 
spite dissipation of vorticity from the numerical code, the vortices remain long-lived because 
baroclinic vorticity production reinforces the vortices and balances the dissipation. Here the 
"sandwich" pattern about each vortex is clearly seen: temperature perturbations track each 
vortex with a warm band to the outside and a cool band to the inside. 



max pert vorticity, rad/yr Kinetic Energy, J/m^ max pert Temperature, K 




0.1 1 10 100 0.1 1 10 100 0.1 1 10 100 

time, orb per time, orb per time, orb per 



Fig. 5. — Comparison of maximum perturbation vorticity \('\ (left), perturbation kinetic 
energy (center), and maximum perturbation temperature \T'\ (right) for simulation Bl, 
where the baroclinic term is turned off at the times indicated. When this occurs the vorticity 
and kinetic energy immediately drop off, indicating that vortex growth is due to the baroclinic 
term. 
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max pert vorticity, rad/yr Kinetic Energy, J/m^ max pert Temperature, K 




0.1 1 10 100 0.1 1 10 100 0.1 1 10 100 

time, orb per time, orb per time, orb per 



Fig. 6. — Data from simulations Taul through Tau4, where r, the radiative coohng time, 
varies between 1 and 100. Two distinct stages can be seen: during vortex formation — at 
early times — the disk cools rapidly, and cools fastest with small r; once vortices have formed 
the baroclinic feedback takes effect, and vortices grow fastest with small r. 
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pert Vort, t=85 orb per pert Vort, t=85 orb per pert Vort, t=85 orb per 




Fig. 7. — Perturbation vorticity and temperature at 85 orbital periods from simulations 
Taul, Tau2, and Tau3, where the radiative cooling time r is varied from 1 to 10. When 
the radiative cooling time is slow (large r, right), the temperature responds sluggishly to 
vortices, making the baroclinic feedback weak. 
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max pert vorticity, rad/yr Kinetic Energy, J/m^ max pert Temperature, K 




1 10 100 1 10 100 0.1 1 10 100 

time, orb per time, orb per time, orb per 



Fig. 8. — Data from simulations Pel-Pe3, which compare the effects of varying Peclet 
number Pe. High Peclet number indicates low thermal diffusion. Increasing thermal diffusion 
(decreasing Pe) strengthens the baroclinic feedback, as exemplified by the slopes of the 
kinetic energy after t = 10. This is similar to increasing the radiative cooling rate. 
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pert Vort, t=50 orb per pert Vort, t=50 orb per 




Fig. 9. — Perturbation vorticity and temperature after 50 orbital periods for simulations 
Pe2 and Pe3. Simulation Pe3 (right) has higher thermal diffusion (lower Pe), resulting in 
temperature perturbations that are larger in size (bottom right), and a stronger baroclinic 
feedback. 



pert Vort, t=60 orb per pert Vort, t=60 orb per pert Vort, t=60 orb per 




Fig. 10. — Perturbation vorticity for simulations Tl, T2, and T3, where the background 
temperature Tq ~ r and d = —2, —1, and —0.75. Larger background temperature gradients 
produce stronger baroclinic instabilities, so that vortices grow faster. 
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max pert vorticity, rad/yr Kinetic Energy, J/m^ max pert Temperature, K 




0.1 1 10 100 0.1 1 10 100 0.1 1 10 100 

time, orb per time, orb per time, orb per 



Fig. 11. — Comparison of data for simulations T1-T5, where the background temperature 
To ~ r'^, and d ranges from -0.25 to -2. Simulations where the background temperature 
gradient is larger in magnitude increases the strength of the baroclinic feedback, resulting in 
increases in all three measures. 



max pert vorticity, rad/yr Kinetic Energy, J/n? max pert Temperature, K 




0.1 1 10 100 0.1 1 10 100 0.1 1 10 100 

time, orb per time, orb per time, orb per 



Fig. 12. — Comparison of data for simulations D1-D3, where the background surface density 
So ~ r ^ and b = —1, —3/2, and —2. This data shows that varying the gradient of surface 
density results in little difference in the strength of the baroclinic feedback. 
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Fig. 13. — Richardson number, Ri{r), for simulations Tl through T5, where the background 
temperature Tq ~ r The Richardson number depends only on background functions, and 
so is constant in time. By comparing to Fig. [H], one can see that a more negative Ri (exp. 
Tl, d = —2) indicates a stronger instability, and less negative Ri (exp. T3, d = —0.75) 
indicates weaker instability. For exp. T4 [d = —0.5) and T5 {d = 0.25) Ri is zero and 
positive, respectively, and the baroclinic instability is not active in both cases. 




5 6 7 8 9 10 



radial distance, au 

Fig. 14. — Azimuthal averages of a/ Ma? for a typical simulation (Al) for various times, in 
orbital periods. This quotient ranges between 10~^ and 10~^ for all simulations. 
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Fig. 15. — Vorticity from the analytic example (top) and numerical model (bottom) where 
vortices are radially leaning out (left) or leaning in (right). The inset shows the angular 
momentum flux, which is positive for outward leaning vortices and negative for inward 
leaning vortices. For these pedagogical examples, the scale for the vorticity and angular 
momentum flux is arbitrary. 
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Fig. 16. — Global angular momentum transport (left) and angular momentum transport 
as a function of radius (right) for simulation Al, where each curve is an average of ten 
measurements taken over one-half orbital period, and times correspond to the images shown 
in Fig. [31 As described in the text, angular momentum transport due to a particular vortex 
depends on the orientation of the vortices. The spatial and temporal variability in angular 
momentum transport shown here is due to variability in the orientation of the vortices. 
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name 


grid 


r 


d 


b 


Re 


Pe 


endtime 


Al 


512^ 


1 


-2 


-3/2 


4e7 


4e7 


100 


A2 


512^ 


1 


-1 


-3/2 


4e7 


4e7 


60 


A3 


5122 


1 


-0.75 


-3/2 


4e7 


4e7 


60 


Bl 


2562 


1 


-0.75 


-3/2 


2e7 


2e7 


300 


Taul 


256^ 


1 


-2 


-3/2 


2e7 


2e7 


300 


Tau2 


2562 


3 


-2 


-3/2 


2e7 


2e7 


300 


Tau3 


256^ 


10 


-2 


-3/2 


2e7 


2e7 


300 


Tau4 


256^ 


100 


-2 


-3/2 


2e7 


2e7 


300 


Tl 


2562 


1 


-2 


-3/2 


2e7 


2e7 


300 


T2 


2562 


1 


-1 


-3/2 


2e7 


2e7 


300 


T3 


2562 


1 


-0.75 


-3/2 


2e7 


2e7 


600 


T4 


2562 


1 


-0.5 


-3/2 


2e7 


2e7 


300 


T5 


2562 


1 


-0.25 


-3/2 


2e7 


2e7 


300 


Dl 


2562 


1 


-2 


-1 


2e7 


2e7 


300 


D2 


2562 


1 


-2 


-3/2 


2e7 


2e7 


300 


D3 


2562 


1 


-2 


-2 


2e7 


2e7 


300 


Pel 


2562 


100 


-2 


-3/2 


2e7 


2e7 


300 


Pe2 


2562 


100 


-2 


-3/2 


2e7 


le6 


300 


Pe3 


2562 


100 


-2 


-3/2 


2e7 


le4 


300 



Table 1: Model parameters for the numerical simulations discussed in this paper. Here r is 
the radiative cooling time in orbital periods, and d and b are the powers on the background 
temperature and surface density functions, Re and Pe are the Reynolds and Peclet numbers, 
and endtime is in orbital periods. In simulation Bl, the baroclinic term is turned off at 
various times during the simulation. 



